clear;

%% Test Ah times a vector
% Set parameters (single tx)

invP = invParams;

%Mesh
nc = 4;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;

% Make a model

k = 1/invP.disT(1,2);

m = abs(rand(invP.mesh.nc,1))/2 + 0.1;

T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - m)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - m)));

M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

% Build Ah

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];

% Test Ah times a vector

k = 1/invP.disT(1,2);
for ii = 1:invP.nt
    j{ii} = rand(invP.mesh.nf,invP.ntx);
end

h = [];
for ii = 1:invP.nt
    h{ii} = rand(invP.mesh.ne,invP.ntx);
end

e = [];
for ii = 1:invP.nt
    e{ii} = rand(invP.mesh.nf,invP.ntx);
end

u = [];
for ii = 1:invP.nt
    u = [u; h{ii}; e{ii}; j{ii}];
end



Ahat_vec1 = Ah*u;
[hq,eq,jq] = calcAhVec(invP,m,h,e,j);
Ahat_vec2 = [cell2mat(hq');cell2mat(eq');cell2mat(jq')]; 
Ahat_vec2 = Ahat_vec2(:);

figure(1);
subplot(2,1,1);
plot(Ahat_vec1,'o');
hold on;
plot(Ahat_vec2,'rx');
hold off;
subplot(2,1,2);
plot(abs(Ahat_vec1-Ahat_vec2)./abs(Ahat_vec1));


%% Test v = A\q
clear
invP = invParams;

%Mesh
nc = 6;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = abs(rand(invP.mesh.nc,1))/2 + 0.1;

T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - m)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - m)));

M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

% Build Ah

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];


q3 = [];
for ii = 1:invP.nt
    q3{ii,1} = rand(invP.mesh.nf,invP.ntx);
end

q = [zeros(invP.mesh.ne+invP.mesh.nf,invP.nt); cell2mat(q3')];
q = q(:);

u_rec1 = Ah\q;
u_rec1 = reshape(u_rec1,invP.mesh.ne+2*invP.mesh.nf,invP.nt);
h_rec1 = u_rec1(1:invP.mesh.ne,:);
h_rec1 = h_rec1(:);

h_rec2 = calcAhInv(invP,m,q3);
h_rec2 = cell2mat(h_rec2);


figure(1);
subplot(2,1,1);
plot(h_rec1,'o');
hold on;
plot(h_rec2,'rx')
hold off
subplot(2,1,2);
plot(abs(h_rec1-h_rec2)./abs(h_rec1));

%%  Test forward modelling

clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = abs(rand(invP.mesh.nc,1))/2 + 0.1;

T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - m)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - m)));

M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

% Build Ah

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];

%set up u0
u0 = [invP.h0; zeros(2*invP.mesh.nf,1)];
q = [-B*u0; zeros((invP.nt-1)*(invP.mesh.ne+2*invP.mesh.nf),1)];

%Solve matrix system

u1 = Ah\q;

[hu2,eu2,ju2] = calcHs_TDIP(invP,m,false);

u2 = [cell2mat(hu2'); cell2mat(eu2'); cell2mat(ju2')];
u2 = u2(:);

figure(1);
subplot(2,1,1);
plot(u1,'o');
hold on;
plot(u2,'rx');
hold off;
subplot(2,1,2);
plot(abs(u1-u2)./abs(u1+1e-8));

%%  Derivative Test derivation of G

clear
invP = invParams;

%Mesh
nc = 6;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = rand(invP.mesh.nc,1);
dm = rand(invP.mesh.nc,1);

% Random fields

j = [];
for ii = 1:invP.nt
    j{ii} = rand(invP.mesh.nf,invP.ntx);
end

h = [];
for ii = 1:invP.nt
    h{ii} = rand(invP.mesh.ne,invP.ntx);
end

e = [];
for ii = 1:invP.nt
    e{ii} = rand(invP.mesh.nf,invP.ntx);
end

u = [];
for ii = 1:invP.nt
    u = [u; h{ii}; e{ii}; j{ii}];
end

fT = @(x)sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - x)));
fTt = @(x)sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - x)));
M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

% Build matricies

T = fT(m);
Tt = fTt(m);

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah_m = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];

T = fT(m+dm);
Tt = fTt(m+dm);

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah_mdm = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];

% Test with matrix system

G = [zeros(invP.mesh.ne,invP.mesh.nc); 
    zeros(invP.mesh.nf,invP.mesh.nc);
    diag(j{1})*invP.mesh.Af*diag(k*invP.tauModel);
    zeros(invP.mesh.ne,invP.mesh.nc); 
    zeros(invP.mesh.nf,invP.mesh.nc);
    diag(j{2} - j{1})*invP.mesh.Af*diag(k*invP.tauModel);
    zeros(invP.mesh.ne,invP.mesh.nc); 
    zeros(invP.mesh.nf,invP.mesh.nc);
    diag(j{3} - j{2})*invP.mesh.Af*diag(k*invP.tauModel)];

norm(Ah_mdm*u - Ah_m*u)
norm(Ah_mdm*u - Ah_m*u - G*dm)

%%  Test G*vec code

clear
invP = invParams;

%Mesh
nc = 6;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = rand(invP.mesh.nc,1);

% Random fields

j = [];
for ii = 1:invP.nt
    j{ii} = rand(invP.mesh.nf,invP.ntx);
end

G = [zeros(invP.mesh.ne,invP.mesh.nc); 
    zeros(invP.mesh.nf,invP.mesh.nc);
    diag(j{1})*invP.mesh.Af*diag(k*invP.tauModel);
    zeros(invP.mesh.ne,invP.mesh.nc); 
    zeros(invP.mesh.nf,invP.mesh.nc);
    diag(j{2} - j{1})*invP.mesh.Af*diag(k*invP.tauModel);
    zeros(invP.mesh.ne,invP.mesh.nc); 
    zeros(invP.mesh.nf,invP.mesh.nc);
    diag(j{3} - j{2})*invP.mesh.Af*diag(k*invP.tauModel)];

Gm1 = G*m;
Gm1 = reshape(Gm1,invP.mesh.ne+2*invP.mesh.nf,invP.nt);
Gm1 = Gm1(end-invP.mesh.nf+1:end,:);
Gm1 = Gm1(:);

Gm2 = Geta_vec(invP,j,m);
Gm2 = cell2mat(Gm2);
Gm2 = Gm2(:);

figure(1);
subplot(2,1,1);
plot(Gm1,'o');
hold on;
plot(Gm2,'rx');
hold off;
subplot(2,1,2);
plot(abs(Gm1-Gm2)./abs(Gm1));

%% Derivative test for dh/dm


clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = rand(invP.mesh.nc,1);
dm = rand(invP.mesh.nc,1);
%%
h = 0.01;

[hm,~,jm] = calcHs_TDIP(invP,m,false);
[hmdm,~,~] = calcHs_TDIP(invP,m + h*dm,false);

Gdm = Geta_vec(invP,jm,dm);
Jdm = calcAhInv(invP,m,Gdm);


hm = cell2mat(hm);
hmdm = cell2mat(hmdm);
Jdm = -cell2mat(Jdm);


norm(hmdm - hm)
norm(hmdm - hm - h*Jdm)


%% Derivative test for J


clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [5 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = rand(invP.mesh.nc,1);
dm = rand(invP.mesh.nc,1);
%%
h = 0.001;

[hm,~,jm] = calcHs_TDIP(invP,m,false);
[hmdm,~,~] = calcHs_TDIP(invP,m + h*dm,false);

Gdm = Geta_vec(invP,jm,dm);
Jdm = calcAhInv(invP,m,Gdm);


dat_m = projHs(invP,hm);
dat_mdm = projHs(invP,hmdm);
J_m = -projHs(invP,Jdm);

norm(dat_mdm - dat_m)
norm(dat_mdm - dat_m - h*J_m)


%% Derivative test for J (function)


clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [5 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = rand(invP.mesh.nc,1);
dm = rand(invP.mesh.nc,1);
%%
h = 0.00001;

[hm,~,jm] = calcHs_TDIP(invP,m,false);
[hmdm,~,~] = calcHs_TDIP(invP,m + h*dm,false);

dat_m = projHs(invP,hm);
dat_mdm = projHs(invP,hmdm);
J_m = Jvec_TDIP(invP,m,jm,dm);

norm(dat_mdm - dat_m)
norm(dat_mdm - dat_m - h*J_m)


%% Adjoint test for Gvec and GTvec codes

clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [5 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;


%Random currents
j = [];
for ii = 1:invP.nt
    j{ii} = rand(invP.mesh.nf,invP.ntx);
end

%random model
w = rand(invP.mesh.nc,1);
v = [];
for ii = 1:invP.nt
    v{ii} = rand(invP.mesh.nf,invP.ntx);
end
vv = cell2mat(v);
vv = vv(:);

t1 = Geta_vec(invP,j,w);
t1 = cell2mat(t1);
t1 = t1(:);
t1 = dot(t1,vv);

t2 = GetaT_vec(invP,j,v);
t2 = dot(t2,w);

%% Test the derivation of hte adjoint problem 

clear
invP = invParams;

%Mesh
nc = 6;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = abs(rand(invP.mesh.nc,1))/2 + 0.1;

T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - m)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - m)));

M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

% Build Ah

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];

qh = [];
for ii = 1:invP.nt
    qh{ii} = rand(invP.mesh.ne,invP.ntx);
end

qe = [];
for ii = 1:invP.nt
    qe{ii} = rand(invP.mesh.nf,invP.ntx);
end

qj = [];
for ii = 1:invP.nt
    qj{ii} = rand(invP.mesh.nf,invP.ntx);
end

q = [];
for ii = 1:invP.nt
    q = [q; qh{ii}; qe{ii}; qj{ii}];
end

u1 = (Ah')\q;
u1 = reshape(u1,invP.mesh.nf*2 + invP.mesh.ne,invP.nt);

h = u1(1:invP.mesh.ne,:);
e = u1(invP.mesh.ne+1:invP.mesh.ne+invP.mesh.nf,:);
j = u1(invP.mesh.ne+invP.mesh.nf+1:end,:);


% figure(1);
% plot(j(:,end),'o');
% hold on;
% plot(-T\e(:,end) - T\qj{end},'rx');
% hold off

% figure(2);
% plot(e(:,end),'o');
% hold on;
% tmp = (S\T)*invP.mesh.CURL*h(:,end) - qj{end} - (S\T)*qe{end};
% plot(tmp,'rx');
% hold off;

% tmp1 = (invP.mesh.CURL'*(S\T)*invP.mesh.CURL + M) * h(:,end);
% tmp2 = qh{end} + invP.mesh.CURL'*(S\T)*qe{end} + invP.mesh.CURL'*qj{end};
% figure(3);
% plot(tmp1,'o');
% hold on;
% plot(tmp2,'rx');
% hold off

% n = 1;
% tmp = -T\e(:,n) + (T\Tt)*j(:,n+1) - T\qj{n};
% figure(4);
% plot(j(:,n),'o');
% hold on;
% plot(tmp,'rx');
% hold off


% n = 2;
% tmp = (S\T)*invP.mesh.CURL*h(:,n) + (Tt - (S\T)*St)*j(:,n+1) - (S\T)*qe{n} - qj{n};
% figure(4);
% plot(e(:,n),'o');
% hold on;
% plot(tmp,'rx');
% hold off


n = 2;
tmp1 = (invP.mesh.CURL'*(S\T)*invP.mesh.CURL + M)*h(:,n);
tmp2 = M*h(:,n+1)+invP.mesh.CURL'*((S\T)*St - Tt)*j(:,n+1) + qh{n} ...
    +invP.mesh.CURL'*(S\T)*qe{n} + invP.mesh.CURL'*qj{n};
figure(4);
plot(tmp1,'o');
hold on;
plot(tmp2,'rx');
hold off


%% Test the adjoint problem 

clear
invP = invParams;

%Mesh
nc = 6;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;
 
% Make a model

k = 1/invP.disT(1,2);

m = abs(rand(invP.mesh.nc,1))/2 + 0.1;

T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - m)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - m)));

M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

% Build Ah

A = [M invP.mesh.CURL' zeros(invP.mesh.ne,invP.mesh.nf);
    invP.mesh.CURL zeros(invP.mesh.nf,invP.mesh.nf) -speye(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) S -T];

B = [-M zeros(size(invP.mesh.CURL')) zeros(invP.mesh.ne,invP.mesh.nf);
    zeros(size(invP.mesh.CURL)) zeros(invP.mesh.nf) zeros(invP.mesh.nf);
    zeros(invP.mesh.nf,invP.mesh.ne) -St Tt];
Ah = [A zeros(size(A)) zeros(size(A)); B A zeros(size(A)); zeros(size(A)) B A];

qh = [];
for ii = 1:invP.nt
    qh{ii} = rand(invP.mesh.ne,invP.ntx);
end

qe = [];
for ii = 1:invP.nt
    qe{ii} = zeros(invP.mesh.nf,invP.ntx);
end

qj = [];
for ii = 1:invP.nt
    qj{ii} = zeros(invP.mesh.nf,invP.ntx);
end

q = [];
for ii = 1:invP.nt
    q = [q; qh{ii}; qe{ii}; qj{ii}];
end

u1 = (Ah')\q;

u1 = reshape(u1,invP.mesh.nf*2 + invP.mesh.ne,invP.nt);

uj1 = u1(invP.mesh.ne+invP.mesh.nf+1:end,:);
uj1 = uj1(:);

uj2 = calcAhTInv(invP,m,qh);
uj2 = cell2mat(uj2);

figure(1);
plot(uj1,'o');
hold on;
plot(uj2,'rx');
hold off;

%% Adjoint test of Qvec and QTvec


clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2;
    sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [3 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Random fields
h = [];
for ii = 1:invP.nt
    h{ii} = rand(invP.mesh.ne,invP.ntx);
end
hv = cell2mat(h);
hv = hv(:);

%Random data
d = rand(invP.nt,invP.ntx);

Qv = projHs(invP,h);
wtQv = dot(d(:),Qv(:));

Qtw = projTHs(invP,d);
Qtw = cell2mat(Qtw);
Qtw = Qtw(:);
vtQtw = dot(hv,Qtw);

%% Adjoint test of Jvec and JTvec


clear
invP = invParams;

%Mesh
nc = 12;
invP.mesh = meshDM(10*ones(nc,1),10*ones(nc,1),10*ones(nc,1));
invP.activeCells = invP.mesh.CC(:,3)<sum(invP.mesh.hz)/2;

%Tx/Rx locations
invP.loc = [sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2 sum(invP.mesh.hx)/2];

%Time stepping
invP.disT = [10 1e-6];
   
%Primary Fields   
invP.h0 = calcH0_VMD(invP);

%Interpolation matrix
invP.Q = getInterpmat(invP);

%Sigma and tau models
invP.sigmaModel = 1e-8*ones(invP.mesh.nc,1);
invP.sigmaModel(invP.activeCells) = 1e-3;

invP.tauModel = zeros(invP.mesh.nc,1);
invP.tauModel(invP.activeCells) = 0.1;

etaModel = rand(invP.mesh.nc,1);

%Random model vector
m = rand(invP.mesh.nc,1);

%Random data
d = rand(invP.nt,invP.ntx);

%Random current
j = [];
for ii = 1:invP.nt
    j{ii} = rand(invP.mesh.nf,invP.ntx);
end

Jm = Jvec_TDIP(invP,etaModel,j,m);
dtJm = dot(d(:),Jm(:));

Jtd = JTvec_TDIP(invP,etaModel,j,d);
mtJtd = dot(m(:),Jtd(:));

dtJm - mtJtd

